Generalized Dissimilarity Modelling (GDM)

Genomic offset predictions

1 Introduction

Resources on Generalized Dissimilarity Modelling (GDM):

• Papers: Ferrier et al. (2007), Mokany et al. (2022) and Fitzpatrick and Keller (2015).

• GDM website: https://mfitzpatrick.al.umces.edu/gdm/

• GDM Github: https://github.com/fitzLab-AL/gdm

From the GDM website: ‘The R package gdm implements Generalized Dissimilarity Modeling (Ferrier et al. 2007) to analyze and map spatial patterns of biodiversity. GDM models biological variation as a function of environment and geography using distance matrices – specifically by relating biological dissimilarity between sites to how much sites differ in their environmental conditions (environmental distance) and how isolated they are from one another (geographical distance). […] GDM also can be used to model other biological levels of organization, notably genetic (Fitzpatrick and Keller 2015) [..] and the approaches for doing so are largely identical to the species-level case with the exception of using a different biological dissimilarity metric depending on the type of response variable.’

2 Formatting data

2.1 Genomic data

2.1.1 Pairwise \(F_{ST}\) matrices

The GDM analysis takes as inputs matrices of pairwise \(F_{ST}\).

To estimate the pairwise \(F_{ST}\), we use the individual-level (i.e. allele counts for each genotype) genomic data with missing data (i.e. no imputation of the missing data), and without minor allele frequencies.

We use the Weir and Cockerham method, with the function pairwise.WCfst of the hierfstat package. We could have used the function gene.dist (with the option WC84), as in Gougherty, Keller, and Fitzpatrick (2021) (and the preprint Gougherty et al. 2020). gene.dist keeps only the lower triangle of the \(F_{ST}\) matrix, while pairwise.WCfst keeps the whole matrix.

# we load the allele counts of each genotype
geno <- read.csv(here("data/DryadRepo/FormattedFilteredGenomicData_AlleleCounts_withoutmaf.csv"), row.names = 1)

geno[geno ==1] <- 12
geno[geno ==2] <- 22
geno[geno ==0] <- 11

geno <- geno %>% 
  t() %>% 
  as.data.frame() %>% 
  dplyr::mutate(pop=substr(row.names(.), 0, 3)) %>% 
  dplyr::select(pop, everything())

geno[1:10,1:10] %>% kable_mydf(boldfirstcolumn = T)

	pop	snp_2	snp_3	snp_4	snp_5	snp_6	snp_7	snp_8	snp_9	snp_10
ALT10	ALT	11	11	11	11	11	11	11	11	11
ALT2	ALT	11	11	11	11	11	12	11	11	11
ALT3	ALT	11	11	11	11	11	11	12	11	11
ALT4	ALT	11	11	12	11	11	12	11	11	11
ALT5	ALT	11	11	11	11	11	11	11	11	11
ALT6	ALT	12	11	12	11	11	11	11	11	11
ALT8	ALT	11	11	11	11	11	12	11	11	11
ALT9	ALT	11	11	11	11	11	11	11	11	11
ARM1	ARM	11	11	11	11	11	11	11	12	11
ARM10	ARM	11	12	11	11	11	11	11	11	11

We load the sets of candidate and control SNPs.

snp_sets <- readRDS(here("outputs/list_snp_sets.rds"))

We calculate the pariwise \(F_{ST}\) matrices.

fst_matrices <- sapply(snp_sets, function(x){ 
  
geno %>% 
  dplyr::select(pop,all_of(x$set_snps)) %>%
  pairwise.WCfst(diploid=TRUE)
    
}, USE.NAMES = TRUE,simplify=FALSE)

# save it
saveRDS(fst_matrices,file=here("data/GDManalysis/FstMatrices.rds"))

2.1.2 Negative values

Pairwise \(F_{ST}\) matrices contains some negative values.

neg_vals <- sapply(fst_matrices, function(x){ length(x[which(x <0)])}, USE.NAMES = TRUE,simplify=FALSE)

lapply(names(snp_sets), function(x){
  
  tibble('SNP set' = snp_sets[[x]]$set_name,
         'Nb of negative values' = neg_vals[[x]])
  
}) %>% 
  bind_rows() %>% 
  kable_mydf()

SNP set	Nb of negative values
All candidate SNPs (380)	6
Common candidate SNPs (69)	34
Candidate SNPs considering pop. struct. correction (221)	8
Control SNPs unmatching allele frequencies (380)	6
Control SNPs matching allele frequencies (380)	8
All SNPs (9817)	2

It generally means there is more variation within than among populations and is likely to result from uneven sample sizes.

2.1.3 Scaling the \(F_{ST}\) matrices

We have to scale the \(F_{ST}\) matrices between 0 and 1 to facilitate model convergence and to enable comparisons the different sets of SNPs, which display different ranges of observed \(F_{ST}\) values.

We standardize the \(F_{ST}\) values in the following way:

\[x_{new} = \frac{x- x_{min}}{x_{max}-x_{min}} \]

This was proposed here: https://www.statisticshowto.com/probability-and-statistics/normal-distributions/normalized-data-normalization/

fst_matrices <- sapply(fst_matrices, function(x){
  
 (x-min(x,na.rm=T))/(max(x,na.rm=T)-min(x,na.rm=T))
  
}, USE.NAMES = TRUE,simplify=FALSE)


# To check that it worked
# sapply(fst_matrices, function(x){range(x,na.rm=T)}, USE.NAMES = TRUE,simplify=FALSE)

2.1.4 Matrix visualization

source(here("scripts/functions/corpmat.R"))
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))

viz_matrix <- function(x,plot_title) {

  p.mat <- corpmat(x)

  p <- corrplot::corrplot(x, method="color", col=col(200),
                   type="upper", order="hclust", 
                   addCoef.col = "black", # Add coefficient of correlation
                   tl.col="black", tl.srt=23, #Text label color and rotation
                   # Combine with significance
                   p.mat = p.mat, sig.level = 0.05, insig = "blank", number.cex =0.8,tl.cex = 0.8,
                   # hide correlation coefficient on the principal diagonal
                   diag=FALSE,
                   title = plot_title,mar=c(0,0,2,0)) # add an upper margin to see the title
  
  return(p)
}

lapply(names(snp_sets), function(x) viz_matrix(x=fst_matrices[[x]], 
                                   plot_title = snp_sets[[x]]$set_name))

2.1.5 Matrix formatting for GDM

For the gdm package: the distance matrix must have as the first column the names of the populations.

fst_matrices <- sapply(fst_matrices, function(x){
  
x <- x %>% 
    as.data.frame() %>% 
    tibble::rownames_to_column("pop") 
  
}, USE.NAMES = TRUE,simplify=FALSE)


saveRDS(fst_matrices,file=here("data/GDManalysis/ScaledFstMatrices.rds"))

fst_matrices <- readRDS(here("data/GDManalysis/ScaledFstMatrices.rds"))

2.2 Climatic data

We do not have to scale the climatic data before the GDM analysis, as scaling of predictors is part of model fitting in GDM, as said in Mokany et al. (2022): ‘As different predictors are measured on different scales (e.g., temperature in degrees, precipitation in mm), they are transformed as part of model fitting, such that the transformed distance between a pair of sites for different predictors can be meaningfully compared and combined.’

# Selected climatic variables
# ===========================
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))

extract_climatedt_metadata(var_clim = clim_var) %>% 
  dplyr::select(label,description,unit) %>% 
  set_colnames(str_to_title(colnames(.))) %>% 
  kable_mydf(font_size = 12)

Label	Description	Unit
bio1	Mean annual temperature	Celsius degrees (°C)
bio3	Isothermality (bio2/bio7) (×100)	index
bio4	Temperature seasonality (standard deviation ×100)	Celsius degrees (°C)
bio12	Annual precipitation	Millimeters (mm)
bio15	Precipitation seasonality (coefficient of variation)	index
SHM	Summer heat moisture index	°C / mm

# Loading point estimate climatic data
adj <- "noADJ"  # not adjusted for elevation
ref_period <- "ref_1901_1950" # reference period 1901-1950

clim_past <- readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]$ref_means %>%
  dplyr::select(pop,longitude,latitude,any_of(clim_var))

# Loading future climatic data extracted from the rasters for each GCM
list_clim_fut <- readRDS(here("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationValuesExtractedFromRasters_FiveGCMs_2041-2070_SSP370.rds")) %>% 
  lapply(function(clim_fut){clim_fut %>% dplyr::select(pop,gcm,any_of(clim_var))})

Warning! (from the GDM website) Note that if your site coordinates are longitude-latitude, the calculation of geographic distances between sites will have errors, the size of which will depend on the geographic extent and location of your study region. We hope to deal with this in a later release, but for now you can avoid these problems by using a projected coordinate system (e.g., equidistant).

Therefore, we reproject the geographic (spheroid) CRS of the population coordinates (with units in degrees longitude and latitude) to a projected (two-dimensional; cartesian) CRS (typically with units of meters from a datum). We will use the CRS EPSG:3035.

# Steps:
  # extract the geographic coordinates of the populations
  # transform the coordinates in spatial points and specify the CRS (WGS84) with sp package
  # transform to a SpatVector object for the terra package
  # reproject with terra package in CRS EPSG:3035
  # extract population coordinates from the SpatVector with terra package

pop_coord_proj <- clim_past %>%
  dplyr::select(contains("ude")) %>% 
  sp::SpatialPoints(proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) %>% 
  terra::vect() %>% 
  terra::project("EPSG:3035") %>% 
  terra::crds() %>% 
  as_tibble() %>% 
  dplyr::rename(long_EPSG_3035=x, lat_EPSG_3035=y)

# merge the population coordinates with the climatic variables
clim_past <- bind_cols(clim_past,pop_coord_proj) %>% dplyr::select(-contains("ude"))
list_clim_fut <- lapply(list_clim_fut, function(clim_fut) bind_cols(clim_fut,pop_coord_proj) %>% dplyr::select(-contains("ude")))

3 GDM fitting

Combining genomic and climatic data with formatsitepair

We use the formatsitepair function to combine the genomic and climatic data into population-pair table format. In our case, we use a pairwise biological dissimilarity matrix (i.e. pairwise \(F_{ST}\)) as response variable, and therefore we have to set bioFormat=3 in the formatsitepair function.

Warning! The rows and columns of the distance matrix have to be in the same order as the rows of the climatic variables. And the distance matrix and the climatic data must not include NAs.

GDM fitting

From the GDM website: ’GDM is a nonlinear extension of permutational matrix regression that uses flexible splines and generalized linear modeling (GLM) to accommodate two types of nonlinearity common in ecological datasets: (1) variation in the rate of compositional turnover (non-stationarity) along environmental gradients, and (2) the curvilinear relationship between biological distance and environmental and geographical distance.

Generalized dissimilarity models are fitted with the gdm function.

Different options are available:

• geo=T to specify that the model should be fit with geographical distance.

• splines an optional vector specifying the the number of I-spline basis functions (the default is three, with larger values producing more complex splines).

• knots an optional vector specifying the locations of “knots” along the splines (defaults to 0 (minimum), 50 (median), and 100 (maximum) quantiles when three I-spline basis functions are used). E

From the GDM website: ‘Even though these option are available, using the default values for these parameters will work fine for most applications. In other words, unless you have a good reason, you should probably use the default settings for splines and knots. The effects (and significance) of altering the number of splines and knot locations has not been systematically explored.’

snp_sets <- lapply(snp_sets, function(x){
  
x$gdm_tab <- gdm_tab <- formatsitepair(bioData=fst_matrices[[x$set_code]], 
                                       bioFormat=3,
                                       XColumn="long_EPSG_3035", 
                                       YColumn="lat_EPSG_3035", 
                                       predData=clim_past, 
                                       siteColumn="pop")  

x$gdm_mod <- gdm_mod  <- gdm(data=gdm_tab, geo=TRUE)

return(x)

})

4 GDM plots

Observed genetic distance vs predicted climatic distance

# Plot predicted climatic distance vs observed genetic distance

plots <- lapply(snp_sets, function(x){
  
overlay_x <- seq( from=min(x$gdm_mod$ecological), to=max(x$gdm_mod$ecological), length=length(x$gdm_mod$ecological) )
overlay_y <- 1 - exp( - overlay_x )
dfline <- tibble(x=overlay_x,y=overlay_y)  

tibble(clim=x$gdm_mod$ecological,
       obs=x$gdm_mod$observed) %>% 
  ggplot(aes(x=clim,y=obs)) + 
  geom_point(color=rgb(0,0,1,0.5)) +
  geom_line(data=dfline,aes(x=x,y=y), color="gray60") +
  xlim(c(0,1.6)) + 
  ylim(c(0,1)) + 
  xlab("Predicted climatic distance") +
  ylab("Observed genetic distance") +
  ggtitle(x$set_name) + 
  theme_bw() +
  theme(title = element_text(size=8))
  
})

big_plot <- plot_grid(plotlist = plots, nrow=2)

ggsave(filename = here("figs/GDM/PredictedClimaticDistancevsObservedGeneticDistancePlots.pdf"), 
       big_plot,
       device="pdf", width=10, height=6)

big_plot

Predicted versus observed genetic distance

plots <- lapply(snp_sets, function(x){
  
tibble(obs=x$gdm_mod$observed, # can also be extracted with x$gdm_tab$distance
       pred=x$gdm_mod$predicted) %>% 
  ggplot(aes(x=obs,y=pred)) + 
  geom_abline(intercept = 0, slope = 1, color="gray60") +
  geom_point(color=rgb(0,0,1,0.5)) +
  xlim(c(0,1)) + ylim(c(0,1)) + 
  xlab("Observed genetic distance") +
  ylab("Predicted genetic distance") +
  ggtitle(x$set_name) + 
  theme_bw()+
  theme(title = element_text(size=8))
  
})

big_plot <- plot_grid(plotlist = plots, nrow=2)

ggsave(filename = here("figs/GDM/PredictedvsObservedGeneticDistancePlots.pdf"), 
       big_plot,
       device="pdf", width=10, height=6)

big_plot

I-splines

From the GDM website: ‘The fitted I-splines provide an indication of how population genetic composition changes along each climatic gradient. They are one of the most informative components of a fitted GDM and so plotting and scrutinizing the splines is a major part of interpreting GDM and the analyzed biological patterns.’

I-splines are shown for predictors with non-zero coefficients, i.e. predictors with a relationship with the genetic distance.

I-spline interpretation: (from the GDM website) ’The maximum height of each spline indicates the magnitude of total genetic change along that gradient and thereby corresponds to the relative importance of that predictor in contributing to allelic turnover while holding all other variables constant (i.e., is a partial climatic distance). The spline’s shape indicates how the rate of genetic change varies with position along that gradient. Thus, the splines provide insight into the total magnitude of genetic change as a function of each gradient and where along each gradient those changes are most pronounced.

lapply(snp_sets, function(x){
  
spline_data <- isplineExtract(x$gdm_mod)

# extract variables with importance different from 0
predictors <- which(apply(spline_data$y,2, sum) != 0) %>% names() 

plots <- lapply(predictors, function(predictor){
 
if(predictor=="Geographic"){
  predictor_name <- "Geography (km)"
  spline_data$x[,predictor] <- spline_data$x[,predictor] / 1000
}  else {
predictor_name <- extract_climatedt_metadata(var_clim = clim_var) %>% 
  filter(label %in% predictor) %>% 
  mutate(predictor_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
  pull(predictor_legend)
}
  
tibble(x=spline_data$x[,predictor],
       y=spline_data$y[,predictor]) %>% 
  ggplot() + 
  geom_line(aes(x=x,y=y),color="blue",linewidth=2) +
  ylim(c(0,1.3)) + 
  ylab("Partial genetic distance") +
  xlab(predictor_name) +
  theme(axis.title.x = element_text(size=9)) +
  theme_bw() 
  
})

# make title
title <- ggdraw() + 
  draw_label(x$set_name,
             fontface = 'bold',
             x = 0,
             hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

# merge plots
p <- plot_grid(plotlist=plots)

# If we want to save the figures without title
# ggsave(here(paste0("figs/GDM/IsplinePlot_",x$set_code,"_noTitle.pdf")), p, width=13,height=9, device="pdf")


# merge plots + title
p <- plot_grid(title, p, ncol = 1,rel_heights = c(0.1, 1))

# save the plots
ggsave(here(paste0("figs/GDM/IsplinePlot_",x$set_code,".pdf")), p, width=13,height=9, device="pdf")

})

5 GDM projections

Warning! Note that the formatsitepair function assumes that the coordinates of the sites are in the same coordinate system as the rasters. At present, no checking is performed to ensure this is the case..

6 GDM predictions

From the GDM website: A GDM model can be used to (i) predict the genetic dissimilarity between population pairs in space or between times using the predict function and (ii) transform the predictor variables from their arbitrary climatic scales to a common biological importance scale using the gdm.transform function.

6.1 Genomic offset of the populations

We want to predict the expected degree of maladaptation in units of the response variable (\(F_{ST}\)) for each population.

6.1.1 At the location of the populations

We aim to predict the genomic offset at the location of the studied populations (predictions based on spatial points). For that, we build a dataset with:

• the genetic distance fixed to O. Indeed, we want the disruption of the current gene-climate relationships under future climates, so we assume that the genetic composition is the same between current and future climate to do this calculation.

• weights are fixed to 1.

# dataframe with past climatic data
pop_tab_pred_clim_past <- clim_past %>% 
  dplyr::select(-pop) %>% 
  dplyr::rename(xCoord=long_EPSG_3035,yCoord=lat_EPSG_3035) %>% 
  set_colnames(paste0("s1.",colnames(.)))

# list of dataframe with future climatic data for each GCM
list_pop_tab_pred <- list_clim_fut %>% 
  lapply(function(clim_fut){
  
clim_fut <-  clim_fut %>% 
  dplyr::select(-pop,-gcm) %>% 
  dplyr::rename(xCoord=long_EPSG_3035,yCoord=lat_EPSG_3035) %>% 
  set_colnames(paste0("s2.",colnames(.)))

bind_cols(pop_tab_pred_clim_past, clim_fut) %>% 
  mutate(distance=0,weights=0) %>% 
  dplyr::select(distance,weights,contains("Coord"), everything())})
  
# calculate the genomic offset for each GCM and each snp set
snp_sets <- lapply(snp_sets, function(x){

x$go <- lapply(list_pop_tab_pred, function(pop_tab_pred){
  
predict(x$gdm_mod,pop_tab_pred)

})
return(x)
})

# Function to plot the relationship between the Euclidean climatic distance and the GDM genetic offset
source(here("scripts/functions/make_eucli_plot.R"))

# Calculate the Euclidean climatic distance
source(here("scripts/functions/generate_scaled_clim_datasets.R"))
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))
clim_dfs <- generate_scaled_clim_datasets(clim_var,clim_ref_adj = FALSE)


# Calculate the Euclidean climatic distance
list_dist_env <- clim_dfs$clim_pred %>% lapply(function(clim_pred){
  
Delta = clim_dfs$clim_ref %>% dplyr::select(any_of(clim_var)) - clim_pred %>% dplyr::select(any_of(clim_var)) 
dist_env = sqrt( rowSums(Delta^2) )

})

# list_dist_env <- list_clim_fut %>% lapply(function(clim_fut){
#   
# Delta = clim_past %>% dplyr::select(-pop,-contains("EPSG")) - clim_fut %>% dplyr::select(-pop,-gcm,-contains("EPSG")) 
# dist_env = sqrt( rowSums(Delta^2) )
# 
# })

# Main gene pools (for the figures)
gps <- readRDS(here("data/GenomicData/MainGenePoolPopulations.rds")) %>% 
  arrange(pop)

par(mfrow=c(1,2))


lapply(snp_sets, function(x) {

lapply(names(list_dist_env), function(gcm){
  
make_eucli_plot(
  X = list_dist_env[[gcm]],
  Y = x$go[[gcm]],
  colors = gps$color_main_gp_pop,
  color_names = gps$main_gp_pop,
  ylab = "GDM genomic offset",
  legend_position="topright",
  plot_title = paste0(x$set_name," - ", gcm))

})
}) 

# We generate scatter plots for the Supplementary Information.
# ============================================================

# Axis limits
# ===========
max_go <- lapply(snp_sets, function(z){
  z$go %>% unlist()
}) %>% unlist() %>% max()

range_eucli <- list_dist_env %>% unlist() %>% range()

# Run the function
# ================
lapply(snp_sets, function(set_i) {

p <- lapply(names(list_dist_env), function(gcm){
  
make_ggscatterplot(
  x = list_dist_env[[gcm]],
  y = set_i$go[[gcm]],
  title = gcm,
  range_eucli = range_eucli,
  max_go = max_go)

})

# remove y-labels to graphs in the second column
p[[2]] <- p[[2]] + ylab("")
p[[4]] <- p[[4]] + ylab("")


# remove x-labels to graphs in the second and third rows
p[[1]] <- p[[1]] + xlab("")
p[[2]] <- p[[2]] + xlab("")
p[[3]] <- p[[3]] + xlab("")

p[[6]] <- get_legend(p[[1]])

for(i in 1:5){p[[i]] <- p[[i]]  +  theme(legend.position = "none")} 


p_grid <- plot_grid(plotlist=p, nrow = 3) 

# We save the figure in pdf
p_grid %>% 
  ggsave(here(paste0("figs/GDM/ScatterPlotEucliDistance_",set_i$set_code,".pdf")), 
         .,
         width=7,
         height=8,
         device="pdf")

# we save the figure in png
p_grid %>% 
  ggsave(here(paste0("figs/GDM/ScatterPlotEucliDistance_",set_i$set_code,".png")), 
         .,
         width=7,
         height=8)

})

$all_cand
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_all_cand.png"

$common_cand
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_common_cand.png"

$corrected_cand
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_corrected_cand.png"

$randomfreq_control
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_randomfreq_control.png"

$simfreq_control
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_simfreq_control.png"

$all
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/GDM/ScatterPlotEucliDistance_all.png"

# Function to make the genomic offset maps
source(here("scripts/functions/make_go_map.R"))

# Population coordinates
pop_coord <-  readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]$ref_means %>%
  dplyr::select(pop,longitude,latitude)


# Generate the maps for each set of SNPs and each GCM
lapply(snp_sets, function(x) {

go_maps <- lapply(names(list_clim_fut), function(gcm){
  
df <- pop_coord %>% mutate(GO = x$go[[gcm]])

make_go_map(df = df,
            plot_title = gcm,
            go_limits = c(0, max(x$go[[gcm]])),
            legend_box_background = "white",
            point_size = 3)

})

legend_maps  <- get_legend(go_maps[[1]])

go_maps <- lapply(go_maps, function(y) y + theme(legend.position = "none"))

go_maps$legend_maps <- legend_maps

go_maps <-plot_grid(plotlist=go_maps)

ggsave(here(paste0("figs/GDM/GOMaps_PopLocations_",x$set_code,".pdf")), go_maps, width=10,height=6, device="pdf")
ggsave(here(paste0("figs/GDM/GOMaps_PopLocations_",x$set_code,".png")), go_maps, width=10,height=6)


# =========
# Add title
# =========
title <- ggdraw() + 
  draw_label(
    x$set_name,
    fontface = 'bold',
    x = 0,
    hjust = 0
  ) +
  theme(plot.margin = margin(0, 0, 0, 7))

# merge title and plots
plot_grid(
  title, go_maps,
  ncol = 1,
  rel_heights = c(0.1, 1))

  })

For each GCM, we attribute the value 1 to the top five populations with the highest genomic offset and we attribute the value 0 to the other populations. We then count the number of 1 for each population, which gives the table and map below:

selected_SNP_set <- "common_cand"

source(here("scripts/functions/make_high_go_pop_maps.R"))

high_go_pops <- make_high_go_pop_maps(pop_coord=pop_coord,
                                      list_go = snp_sets[[selected_SNP_set]]$go,
                                      ggtitle = paste0("GDM - ",snp_sets[[selected_SNP_set]]$set_name),
                                      nb_id_pop = 5) # number of selected populations

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2 3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

saveRDS(high_go_pops, file = here("outputs/GDM/high_go_pops.rds"))

high_go_pops[[1]] %>% kable_mydf

pop	longitude	latitude	GFDL-ESM4	IPSL-CM6A-LR	MPI-ESM1-2-HR	MRI-ESM2-0	UKESM1-0-LL	sum_go
VAL	-4.31	40.52	0	0	0	0	0	0
ARN	-5.12	40.19	0	0	0	0	0	0
CEN	-4.49	40.28	0	0	0	0	0	0
BON	-1.66	39.99	0	0	0	0	0	0
ORI	-2.35	37.53	0	0	0	0	0	0
QUA	-0.36	38.97	0	0	0	0	0	0
OLB	-0.62	40.17	0	0	0	0	0	0
PIA	9.46	42.02	0	0	0	0	0	0
VER	-1.09	45.55	0	0	0	0	0	0
MAD	-5.23	35.18	0	0	0	0	0	0
COC	-4.50	41.25	0	0	0	0	0	0
CUE	-4.48	41.37	0	0	0	0	0	0
COM	-3.95	36.83	0	0	0	0	0	0
TAM	-5.02	33.60	0	0	0	0	0	0
HOU	-1.15	45.18	0	0	0	0	0	0
BAY	-2.88	41.52	0	0	0	0	0	0
CAR	-4.28	41.17	0	0	0	0	0	0
STJ	-2.03	46.76	0	0	0	0	0	0
PIE	9.04	41.97	0	0	0	0	0	0
OLO	-1.83	46.57	0	0	0	0	0	0
MIM	-1.30	44.13	0	0	0	0	0	0
PLE	-2.34	47.78	0	0	0	0	0	0
SAL	-3.06	41.83	0	0	0	0	0	0
SIE	-6.49	43.53	0	0	0	0	1	1
LEI	-8.96	39.78	0	0	0	1	0	1
PET	-1.30	44.06	0	1	0	0	0	1
CAS	-6.98	43.50	0	0	0	0	1	1
PUE	-6.63	43.55	0	0	0	0	1	1
LAM	-6.22	43.56	0	1	1	0	0	2
SEG	-8.45	42.82	1	0	0	1	0	2
SAC	-8.36	42.12	1	0	1	1	0	3
CAD	-6.42	43.54	1	1	1	0	1	4
ALT	-6.49	43.28	1	1	1	1	0	4
ARM	-6.46	43.30	1	1	1	1	1	5

high_go_pops[[2]]

6.1.1.1 Comparing GO predictions

We look at the correlation across the different genomic offset predictions at the location of the populations, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

lapply(names(snp_sets[[1]]$go), function(gcm){
  
lapply(snp_sets, function(x) x$go[[gcm]]) %>% 
  as_tibble() %>%
  cor() %>% 
  corrplot(method = 'number',type = 'lower', 
           diag = FALSE,mar=c(0,0,2,0),
           title=gcm,
           number.cex=2,tl.cex=1.5)
  
})

6.1.2 Across the species range

To predict the genomic offset across the species range, we use climatic rasters for the reference period (1901-1950). These rasters were provided by Maurizio Marchi and come from the Climate Downscaling Tool (ClimateDT). The values of climatic variables at the location of the populations match between the rasters and the spatial points (I checked it in the report 3_CheckingPastFutureClimatesPopulationLocations).

6.1.2.1 Maps with raster package

# Past rasters
path <- here("data/ClimaticData/ClimateDTRasters/1km_1901-1950_Extent-JulietteA/")
tif_paths <- lapply(clim_var, function(x) paste0(path,"/",x,".tif"))
past_rasts <- raster::stack(tif_paths) %>% projectRaster(crs = "EPSG:3035")


list_pred_rasts <- lapply(snp_sets, function(x){ # for each snp set
lapply(names(list_clim_fut), function(gcm){ # for each GCM
  
# Future rasters  
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))
tif_paths <- lapply(clim_var, function(x) paste0(path,"/",x,".tif"))
fut_rasts <- raster::stack(tif_paths)
fut_rasts <- projectRaster(fut_rasts, crs = "EPSG:3035")
  

predict(x$gdm_mod,
        past_rasts, # raster with current climate at 2.5 minutes resolution
        time=TRUE,
        fut_rasts)  # Rasters with future climates at 2.5 minutes resolution

}) %>% setNames(names(list_clim_fut))
})


list_pred_rasts %>% saveRDS(here("outputs/GDM/go_pred_rasters.rds"))

list_pred_rasts <- readRDS(here("outputs/GDM/go_pred_rasters.rds"))

par(mfrow=c(1,2))

names(snp_sets) %>% lapply(function(x){
  
lapply(names(list_clim_fut), function(gcm){

raster::plot(list_pred_rasts[[x]][[gcm]], 
             col=rgb.tables(1000),
             axes=FALSE, box=FALSE,
             main=paste0(snp_sets[[x]]$set_name," - ", gcm))
})
  
})

6.1.2.2 Maps with ggplot2

We project the genomic offset predictions for the set of candidate SNPs and the mean genomic offset across the five GCMs.

MANUSCRIPT FIGURE: This figure corresponds to Figure 6a in the main manuscript.

# Buffer for the maps (maritime pine distribution)
range_buffer = shapefile(here('data/Mapping/PinpinDistriEUforgen_NFIplotsBuffer10km.shp'))

# We project the genomic offset only for the candidate SNPs
snp_set <- snp_sets[["all_cand"]]

# We load the rasters with the climates of the reference period
path <- here("data/ClimaticData/ClimateDTRasters/1km_1901-1950_Extent-JulietteA/")
ref_rasts <- lapply(clim_var, function(x) paste0(path,"/",x,".tif")) %>% 
  raster::stack() %>%
  mask(range_buffer) %>% 
  projectRaster(crs = "EPSG:3035")


df <- lapply(names(list_clim_fut), function(gcm){ # for each GCM
  
# Future rasters  
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))
fut_rasts <- lapply(clim_var, function(x) paste0(path,"/",x,".tif")) %>% 
  raster::stack() %>% 
  mask(range_buffer) %>% 
  projectRaster(crs = "EPSG:3035")
  

predict(snp_set$gdm_mod,
              ref_rasts, # rasters with reference climate at 2.5 minutes resolution
              time=TRUE,
              fut_rasts) %>%  # rasters with future climates at 2.5 minutes resolutio
  projectRaster(crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") %>% 
  rasterToPoints() %>% 
  set_colnames(c("x","y","GO")) %>% 
  as_tibble() %>% 
  mutate(gcm=gcm)

}) %>% 
  bind_rows %>% 
  pivot_wider(values_from = GO, names_from = gcm) %>% 
  mutate(mean_GO = rowMeans(dplyr::select(.,-c(x:y))))


df %>% saveRDS(here("outputs/GDM/go_pred_rasters_2_allcand.rds"))

# Map options
# ===========
point_size = 2
x_limits = c(-10, 15)
y_limits = c(31, 52)

# Country borders
world <- ne_countries(scale = "medium", returnclass = "sf")

# Load the mean GO projections
df <- readRDS(here("outputs/GDM/go_pred_rasters_2_allcand.rds"))

p <- ggplot(data=df) + 
  geom_sf(data = world, fill="gray98") + 
  scale_x_continuous(limits = x_limits) +
  scale_y_continuous(limits = y_limits) + 
  geom_raster(aes(x = x, y = y, fill = mean_GO), alpha = 1) + 
  scale_fill_gradient2(low="blue", mid= "yellow", high="red",
                       midpoint=(max(df$mean_GO)-min(df$mean_GO))/2,
                       limits=c(min(df$mean_GO),max(df$mean_GO)),
                       name = "Genomic offset from GDM") +
  xlab("Longitude") + ylab("Latitude") +
  #ggtitle("Genomic offset predictions with GDM") +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), 
        plot.background = element_blank(), 
        panel.background = element_blank(), 
        legend.position = c(0.8,0.15),
        legend.box.background = element_rect(colour = "gray80"),
        legend.title = element_text(size=10),
        strip.text = element_text(size=11))

p %>% ggsave(here("figs/GDM/GOmeanProjections_AllCandSNPs.pdf"),., width=6,height=6, device="pdf")
p %>% ggsave(here("figs/GDM/GOmeanProjections_AllCandSNPs.png"),., width=6,height=6)
p

6.1.3 Corr btw predictions with rasters or spatial points

We check that the genomic offset predictions at the locations of the populations are correlated between those obtained with the spatial points and those obtained with the rasters.

# checking correlations
names(snp_sets) %>% lapply(function(x){
  
cor_go <- lapply(names(list_clim_fut), function(gcm){
  
go_rast <- raster::extract(list_pred_rasts[[x]][[gcm]], clim_past[,c("long_EPSG_3035","lat_EPSG_3035")])

cor(snp_sets[[x]]$go[[gcm]],go_rast) 
  
}) %>% unlist() 
  
tibble(gcm=names(list_clim_fut),
            cor_go=cor_go)
  
}) %>% 
  setNames(names(snp_sets)) %>% 
  bind_rows(.id="snp_set") %>% 
  kable_mydf()

snp_set	gcm	cor_go
all_cand	GFDL-ESM4	0.99
all_cand	IPSL-CM6A-LR	0.99
all_cand	MPI-ESM1-2-HR	0.99
all_cand	MRI-ESM2-0	0.99
all_cand	UKESM1-0-LL	0.99
common_cand	GFDL-ESM4	0.99
common_cand	IPSL-CM6A-LR	0.99
common_cand	MPI-ESM1-2-HR	0.99
common_cand	MRI-ESM2-0	0.99
common_cand	UKESM1-0-LL	0.99
corrected_cand	GFDL-ESM4	0.99
corrected_cand	IPSL-CM6A-LR	0.99
corrected_cand	MPI-ESM1-2-HR	0.99
corrected_cand	MRI-ESM2-0	0.99
corrected_cand	UKESM1-0-LL	0.99
randomfreq_control	GFDL-ESM4	1.00
randomfreq_control	IPSL-CM6A-LR	1.00
randomfreq_control	MPI-ESM1-2-HR	1.00
randomfreq_control	MRI-ESM2-0	1.00
randomfreq_control	UKESM1-0-LL	1.00
simfreq_control	GFDL-ESM4	0.99
simfreq_control	IPSL-CM6A-LR	0.99
simfreq_control	MPI-ESM1-2-HR	0.99
simfreq_control	MRI-ESM2-0	1.00
simfreq_control	UKESM1-0-LL	0.99
all	GFDL-ESM4	0.99
all	IPSL-CM6A-LR	0.99
all	MPI-ESM1-2-HR	1.00
all	MRI-ESM2-0	1.00
all	UKESM1-0-LL	1.00

This is ok!

7 Validation - NFI plots

Climatic data and coordinates

# Load the climatic data of the NFI plots.
nfi_clim <- readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds")) 

# Reproject the plotcode coordinates in the projection system `EPSG:3035`.
plotcoord_EPSG_3035 <- nfi_clim[[1]] %>%
  dplyr::select(contains("ude")) %>% 
  sp::SpatialPoints(proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) %>% 
  terra::vect() %>% 
  terra::project("EPSG:3035") %>% 
  terra::crds() %>% 
  as_tibble() %>% 
  dplyr::rename(xCoord=x, yCoord=y)

# Keep only the climatic variables of interest and the plotcode coordinates
nfi_clim <- nfi_clim %>% 
  lapply(function(x){
    
x %>% bind_cols(plotcoord_EPSG_3035) %>% 
      dplyr::select(contains("Coord"),any_of(clim_var))})

# Format as GDM inputs
nfi_clim[[1]] <- nfi_clim[[1]] %>% set_colnames(paste0("s1.",colnames(.)))
nfi_clim[[2]] <- nfi_clim[[2]] %>% set_colnames(paste0("s2.",colnames(.)))

nfi_clim <- bind_cols(nfi_clim) %>% 
  mutate(distance=0,weights=0) %>% 
  dplyr::select(distance,weights,contains("Coord"), everything())

Calculate the genomic offset for each snp set

snp_sets <- lapply(snp_sets, function(x){

x$go_nfi <- predict(x$gdm_mod,nfi_clim)

return(x)
})

Map genomic offset predictions at the location of the NFI plots

These figures are included in the Supplementary Information.

source(here("scripts/functions/make_go_map.R"))

lapply(snp_sets, function(x) {

df <- readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))[[1]] %>% 
dplyr::select(contains("ude")) %>% mutate(GO = x$go_nfi) 

p <- make_go_map(df = df, 
                 type="NFI",
                 point_size = 0.5,
                 plot_title =  x$set_name,
                 go_limits = c(0,max(x$go_nfi)),
                 legend_position = c(0.85,0.2), 
                 y_limits = c(35, 51))
  
  
ggsave(here(paste0("figs/GDM/NFI_GOmap_",x$set_code,".pdf")), p, width=6,height=6, device="pdf")
ggsave(here(paste0("figs/GDM/NFI_GOmap_",x$set_code,".png")), p, width=6,height=6)
# p <- p + theme(plot.title = element_blank()) # to remove the title

# Show maps in the Quarto document
# ================================
p  
  })

We look at the correlation across the different genomic offset predictions in the NFI plots, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

lapply(snp_sets, function(x) x$go_nfi) %>% 
  as_tibble() %>%
  cor() %>% 
  corrplot(method = 'number',type = 'lower', diag = FALSE,mar=c(0,0,2,0),
               number.cex=2,tl.cex=1.5)

8 Validation - Common gardens

Climatic data and coordinates

# format climatic data
cg_clim <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>%  
  dplyr::select(any_of(clim_var)) %>% 
  set_colnames(paste0("s2.",colnames(.))) %>%  
  replicate(nrow(pop_tab_pred_clim_past),.,simplify=F) %>% 
  bind_rows() %>% 
  bind_cols(slice(pop_tab_pred_clim_past, rep(1:n(), each = 5))) %>% 
  mutate(s2.xCoord=s1.xCoord, # same coordinates, so that only climatic differences are considered
         s2.yCoord=s1.yCoord,
         distance=0,
         weights=0) %>% 
  dplyr::select(distance, weights, s1.xCoord, s1.yCoord, s2.xCoord, s2.yCoord,contains("s1"), contains("s2"))

Calculate the genomic offset for each snp set

go_df <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% 
  dplyr::select(cg) %>% 
  replicate(nrow(pop_tab_pred_clim_past),.,simplify=F) %>% 
  bind_rows() %>% 
  bind_cols(slice(clim_past[,"pop"],rep(1:n(), each = 5)))

snp_sets <- lapply(snp_sets, function(x){

x$go_cg <- go_df %>% 
  mutate(go=predict(x$gdm_mod,cg_clim)) %>% 
  pivot_wider(names_from = cg, values_from = go)

return(x)
})

Map genomic offset predictions at the locations of the populations

cg_coord <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% dplyr::select(cg,contains("ude"))
cg_names <- unique(cg_coord$cg)

go_maps_cg <- lapply(cg_names, function(cg_name){

p <- lapply(snp_sets, function(x) {

df <- pop_coord %>%
        left_join(x$go_cg[,c("pop",cg_name)], by="pop") %>% 
        dplyr::rename(GO=all_of(cg_name))
  
 p <- make_go_map(df = df,
                  plot_title = paste0(str_to_title(cg_name), " - ",x$set_name),
                  point_size = 3,
                  type="CG",
                  go_limits = c(0,max(x$go_cg[[cg_name]])),
                  cg_coord = filter(cg_coord, cg == cg_name))
 
 
ggsave(filename = here(paste0("figs/GDM/GOmap_",x$set_code,"_",cg_name,".pdf")), device = "pdf",width=6,height=6)
ggsave(filename = here(paste0("figs/GDM/GOmap_",x$set_code,"_",cg_name,".png")), width=6,height=6)

# p + theme(plot.title = element_blank(), legend.position = "none") # to rm the title
 
p
 
  })

plot_grid(plotlist=p,nrow=2)
  
}) %>% setNames(cg_names)

pdf(here("figs/GDM/GOmaps_CGs.pdf"), width=15,height=6)
lapply(go_maps_cg, function(x) x)
dev.off()

# show maps
lapply(go_maps_cg, function(x) x)

We look at the correlation across the different genomic offset predictions in the common gardens, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

lapply(cg_names, function(cg_name){

lapply(snp_sets, function(x) x$go_cg) %>% 

lapply(function(set){
  set[[cg_name]]
}) %>% 
    as_tibble() %>% 
    cor() %>% 
    corrplot(method = 'number',type = 'lower', 
             diag = FALSE,mar=c(0,0,2,0),
             title=str_to_title(cg_name),
             number.cex=2,tl.cex=1.5)

})

9 Saving GO predictions

Let’s save the genomic offset predictions for comparison with the other methods.

snp_sets %>% saveRDS(file=here("outputs/GDM/go_predictions.rds"))

10 Session information

devtools::session_info()

─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2025-01-16
 pandoc   3.1.11 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package           * version   date (UTC) lib source
 cachem              1.0.8     2023-05-01 [1] CRAN (R 4.4.0)
 cellranger          1.1.0     2016-07-27 [1] CRAN (R 4.4.0)
 class               7.3-23    2025-01-01 [4] CRAN (R 4.4.2)
 classInt            0.4-10    2023-09-05 [1] CRAN (R 4.4.0)
 cli                 3.6.2     2023-12-11 [1] CRAN (R 4.4.0)
 cluster             2.1.8     2024-12-11 [4] CRAN (R 4.4.2)
 codetools           0.2-19    2023-02-01 [4] CRAN (R 4.2.2)
 colorRamps        * 2.3.4     2024-03-07 [1] CRAN (R 4.4.0)
 colorspace          2.1-0     2023-01-23 [1] CRAN (R 4.4.0)
 corrplot          * 0.92      2021-11-18 [1] CRAN (R 4.4.0)
 cowplot           * 1.1.3     2024-01-22 [1] CRAN (R 4.4.0)
 DBI                 1.2.2     2024-02-16 [1] CRAN (R 4.4.0)
 devtools            2.4.5     2022-10-11 [1] CRAN (R 4.4.0)
 digest              0.6.35    2024-03-11 [1] CRAN (R 4.4.0)
 doParallel          1.0.17    2022-02-07 [1] CRAN (R 4.4.0)
 dplyr             * 1.1.4     2023-11-17 [1] CRAN (R 4.4.0)
 e1071               1.7-14    2023-12-06 [1] CRAN (R 4.4.0)
 ellipsis            0.3.2     2021-04-29 [1] CRAN (R 4.4.0)
 evaluate            0.23      2023-11-01 [1] CRAN (R 4.4.0)
 fansi               1.0.6     2023-12-08 [1] CRAN (R 4.4.0)
 farver              2.1.2     2024-05-13 [1] CRAN (R 4.4.0)
 fastmap             1.1.1     2023-02-24 [1] CRAN (R 4.4.0)
 forcats           * 1.0.0     2023-01-29 [1] CRAN (R 4.4.0)
 foreach             1.5.2     2022-02-02 [1] CRAN (R 4.4.0)
 fs                  1.6.4     2024-04-25 [1] CRAN (R 4.4.0)
 gdm               * 1.5.0-9.1 2022-12-01 [1] CRAN (R 4.4.0)
 generics            0.1.3     2022-07-05 [1] CRAN (R 4.4.0)
 ggplot2           * 3.5.1     2024-04-23 [1] CRAN (R 4.4.0)
 glue                1.7.0     2024-01-09 [1] CRAN (R 4.4.0)
 gtable              0.3.5     2024-04-22 [1] CRAN (R 4.4.0)
 here              * 1.0.1     2020-12-13 [1] CRAN (R 4.4.0)
 hierfstat         * 0.5-11    2022-05-05 [1] CRAN (R 4.4.0)
 highr               0.10      2022-12-22 [1] CRAN (R 4.4.0)
 hms                 1.1.3     2023-03-21 [1] CRAN (R 4.4.0)
 htmltools           0.5.8.1   2024-04-04 [1] CRAN (R 4.4.0)
 htmlwidgets         1.6.4     2023-12-06 [1] CRAN (R 4.4.0)
 httpuv              1.6.15    2024-03-26 [1] CRAN (R 4.4.0)
 httr                1.4.7     2023-08-15 [1] CRAN (R 4.4.0)
 iterators           1.0.14    2022-02-05 [1] CRAN (R 4.4.0)
 jsonlite            1.8.8     2023-12-04 [1] CRAN (R 4.4.0)
 kableExtra        * 1.4.0     2024-01-24 [1] CRAN (R 4.4.0)
 KernSmooth          2.23-26   2025-01-01 [4] CRAN (R 4.4.2)
 knitr             * 1.46      2024-04-06 [1] CRAN (R 4.4.0)
 labeling            0.4.3     2023-08-29 [1] CRAN (R 4.4.0)
 later               1.3.2     2023-12-06 [1] CRAN (R 4.4.0)
 latex2exp         * 0.9.6     2022-11-28 [1] CRAN (R 4.4.0)
 lattice             0.22-5    2023-10-24 [4] CRAN (R 4.3.1)
 lifecycle           1.0.4     2023-11-07 [1] CRAN (R 4.4.0)
 lubridate         * 1.9.3     2023-09-27 [1] CRAN (R 4.4.0)
 magrittr          * 2.0.3     2022-03-30 [1] CRAN (R 4.4.0)
 MASS                7.3-64    2025-01-04 [4] CRAN (R 4.4.2)
 Matrix              1.7-1     2024-10-18 [4] CRAN (R 4.4.1)
 memoise             2.0.1     2021-11-26 [1] CRAN (R 4.4.0)
 mgcv                1.9-1     2023-12-21 [4] CRAN (R 4.3.2)
 mime                0.12      2021-09-28 [1] CRAN (R 4.4.0)
 miniUI              0.1.1.1   2018-05-18 [1] CRAN (R 4.4.0)
 munsell             0.5.1     2024-04-01 [1] CRAN (R 4.4.0)
 nlme                3.1-166   2024-08-14 [4] CRAN (R 4.4.2)
 pbapply             1.7-2     2023-06-27 [1] CRAN (R 4.4.0)
 permute             0.9-7     2022-01-27 [1] CRAN (R 4.4.0)
 pillar              1.9.0     2023-03-22 [1] CRAN (R 4.4.0)
 pkgbuild            1.4.4     2024-03-17 [1] CRAN (R 4.4.0)
 pkgconfig           2.0.3     2019-09-22 [1] CRAN (R 4.4.0)
 pkgload             1.3.4     2024-01-16 [1] CRAN (R 4.4.0)
 plyr                1.8.9     2023-10-02 [1] CRAN (R 4.4.0)
 profvis             0.3.8     2023-05-02 [1] CRAN (R 4.4.0)
 promises            1.3.0     2024-04-05 [1] CRAN (R 4.4.0)
 proxy               0.4-27    2022-06-09 [1] CRAN (R 4.4.0)
 purrr             * 1.0.2     2023-08-10 [1] CRAN (R 4.4.0)
 R6                  2.5.1     2021-08-19 [1] CRAN (R 4.4.0)
 ragg                1.3.0     2024-03-13 [1] CRAN (R 4.4.0)
 raster            * 3.6-26    2023-10-14 [1] CRAN (R 4.4.0)
 Rcpp                1.0.12    2024-01-09 [1] CRAN (R 4.4.0)
 readr             * 2.1.5     2024-01-10 [1] CRAN (R 4.4.0)
 readxl            * 1.4.3     2023-07-06 [1] CRAN (R 4.4.0)
 remotes             2.5.0     2024-03-17 [1] CRAN (R 4.4.0)
 reshape2            1.4.4     2020-04-09 [1] CRAN (R 4.4.0)
 rlang               1.1.4     2024-06-04 [1] CRAN (R 4.4.0)
 rmarkdown           2.26      2024-03-05 [1] CRAN (R 4.4.0)
 rnaturalearth     * 1.0.1     2023-12-15 [1] CRAN (R 4.4.0)
 rnaturalearthdata * 1.0.0     2024-02-09 [1] CRAN (R 4.4.0)
 rprojroot           2.0.4     2023-11-05 [1] CRAN (R 4.4.0)
 rstudioapi          0.16.0    2024-03-24 [1] CRAN (R 4.4.0)
 scales              1.3.0     2023-11-28 [1] CRAN (R 4.4.0)
 sessioninfo         1.2.2     2021-12-06 [1] CRAN (R 4.4.0)
 sf                * 1.0-16    2024-03-24 [1] CRAN (R 4.4.0)
 shiny               1.8.1.1   2024-04-02 [1] CRAN (R 4.4.0)
 sp                * 2.1-4     2024-04-30 [1] CRAN (R 4.4.0)
 stringi             1.8.4     2024-05-06 [1] CRAN (R 4.4.0)
 stringr           * 1.5.1     2023-11-14 [1] CRAN (R 4.4.0)
 svglite             2.1.3     2023-12-08 [1] CRAN (R 4.4.0)
 systemfonts         1.0.6     2024-03-07 [1] CRAN (R 4.4.0)
 terra             * 1.7-71    2024-01-31 [1] CRAN (R 4.4.0)
 textshaping         0.3.7     2023-10-09 [1] CRAN (R 4.4.0)
 tibble            * 3.2.1     2023-03-20 [1] CRAN (R 4.4.0)
 tidyr             * 1.3.1     2024-01-24 [1] CRAN (R 4.4.0)
 tidyselect          1.2.1     2024-03-11 [1] CRAN (R 4.4.0)
 tidyverse         * 2.0.0     2023-02-22 [1] CRAN (R 4.4.0)
 timechange          0.3.0     2024-01-18 [1] CRAN (R 4.4.0)
 tzdb                0.4.0     2023-05-12 [1] CRAN (R 4.4.0)
 units               0.8-5     2023-11-28 [1] CRAN (R 4.4.0)
 urlchecker          1.0.1     2021-11-30 [1] CRAN (R 4.4.0)
 usethis             2.2.3     2024-02-19 [1] CRAN (R 4.4.0)
 utf8                1.2.4     2023-10-22 [1] CRAN (R 4.4.0)
 vctrs               0.6.5     2023-12-01 [1] CRAN (R 4.4.0)
 vegan               2.6-4     2022-10-11 [1] CRAN (R 4.4.0)
 viridisLite         0.4.2     2023-05-02 [1] CRAN (R 4.4.0)
 withr               3.0.0     2024-01-16 [1] CRAN (R 4.4.0)
 xfun                0.43      2024-03-25 [1] CRAN (R 4.4.0)
 xml2                1.3.6     2023-12-04 [1] CRAN (R 4.4.0)
 xtable            * 1.8-4     2019-04-21 [1] CRAN (R 4.4.0)
 yaml                2.3.8     2023-12-11 [1] CRAN (R 4.4.0)

 [1] /home/juliette/R/x86_64-pc-linux-gnu-library/4.4
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

References

Ferrier, Simon, Glenn Manion, Jane Elith, and Karen Richardson. 2007. “Using Generalized Dissimilarity Modelling to Analyse and Predict Patterns of Beta Diversity in Regional Biodiversity Assessment.” Diversity and Distributions 13 (3): 252–64. https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1472-4642.2007.00341.x.

Fitzpatrick, Matthew C, and Stephen R Keller. 2015. “Ecological Genomics Meets Community-Level Modelling of Biodiversity: Mapping the Genomic Landscape of Current and Future Environmental Adaptation.” Ecology Letters 18 (1): 1–16. https://onlinelibrary.wiley.com/doi/abs/10.1111/ele.12376.

Gougherty, Andrew V, Stephen R Keller, and Matthew C Fitzpatrick. 2021. “Maladaptation, Migration and Extirpation Fuel Climate Change Risk in a Forest Tree Species.” Nature Climate Change 11 (2): 166–71.

Mokany, Karel, Chris Ware, Skipton NC Woolley, Simon Ferrier, and Matthew C Fitzpatrick. 2022. “A Working Guide to Harnessing Generalized Dissimilarity Modelling for Biodiversity Analysis and Conservation Assessment.” Global Ecology and Biogeography 31 (4): 802–21. https://onlinelibrary.wiley.com/doi/full/10.1111/geb.13459.